
\chapter{数值误差}

\section{数字的表达}
我们知道，有理数有无限多个，无理数也有无限多个。
但是计算机能表达的数，只可能是有限多个。
我们人有十个指头，所以发展了 10 进制（如果用16进制，3岁小孩可能抱怨数学太难了），我们表达正整数都用这个形式：
\begin{equation}
	a_0 10^0 + a_1 10^1 + \cdots + a_n 10^n
\end{equation}
而计算机用 2 进制，所以机器们这样表达数字（如果电冰箱和电视机交谈的话）
\begin{equation}
	a_0 2^0 + a_1 2^1 + \cdots + a_n 2^n
\end{equation}
位数是有限的，比如 c++ 中的 int 占 2 个 byte，即 16 个字节。

对于浮点数，我们以单精度浮点数为例，它有 4 个 byte，即 32 个 bit,以如下形式表达数字
\begin{equation}
	(-1)^s 0. a_{-1} a_{-2} a_{-3} ... a_{-23} 2^{exponent},
\end{equation}
其中，s 占一个 bit，exponent 占 8 个 bit，所以数值部分占 23 个bit（一共32位）。
所以，上面的式子可以表示一系列离散的数。
如果你将值 x 赋给变量 a，x 无法用上面的形式严格表出，那么计算机会在自己能表达的数中，取一个离 x 最近的数，我们将它记作 $fl(x)$（$fl$ 是 float 的缩写），赋值给 a，所以实际上发生的是：
\begin{equation}
	a = fl (x),
\end{equation}
那么，$fl (x)$ 相对 x 的相对误差的绝对值为
\begin{equation}
	| \frac{ fl (x) - x }{ fl (x) } | \approx 2^{-23} \approx 10^{-7}
\end{equation}
类似地，双精度数的相对误差约为 $10^{-16}$。

这就叫做 {\bf round off error}，我国人民将它翻译为 舍入误差。

\section{舍入误差：两个相近的数相减}
我们说了，$x$ 在计算机中实际上是 $fl(x)$，那么，如果我们要计算两个双精度数的差 $a = b - c$，则实际上发生的是
\begin{equation}
	fl(a) = fl(b) - fl(c),
\end{equation}
它并不是真的 $a$ 的值，我们自然会关心，误差有多大。
相对误差为
\begin{equation}
	\frac{fl(a)}{a} -1 = \frac{ fl(b) - fl(c) }{ a} -1 = \frac{ b (1+ \epsilon_b) - c(1+\epsilon_c) }{a}
	-1= \frac{ b \epsilon_b - c \epsilon_c }{a},
\end{equation}
其中，$\epsilon_b, \epsilon_c$ 是 $fl(b), fl(c)$ 的相对误差，都是 $10^{-16}$ 量级
如果 $ b \approx c$，我们可以将上式近似为
\begin{equation}
	\frac{ fl(a) }{a} -1= \frac{ b \epsilon_b - c \epsilon_c }{a} \approx \frac{c}{a}(\epsilon_b - \epsilon_c),
\end{equation}
$\epsilon_b, \epsilon_c$ 都是 $2^{-23}$ 数量级，所以它俩的差也约为 $10^{-16}$。
所以对于双精度数，如果 $ a= b-c << c$，相对误差的上限约为
\begin{equation}
	\frac{ fl(a) }{a} -1 \leq \frac{c}{a} 10^{-16}
\end{equation}
相应地，误差的上限约为
\begin{equation}
	fl(a) - a \leq c 10^{-16}
\end{equation}
这个误差有可能超过 $a$，就完全失去了精度，计算失去了意义。

如果用单精度浮点数，上面的 $10^{-16}$ 就变成 $10^{-7}$，所以会更加糟糕。
所以在本讲义中，我们一律用 double 而不用 float。

我们做一个超级简单的例子，让 10000.1 减去 10000，然后命令计算机，输出结果的前 15 位有效数字。
\begin{lstlisting}
#include<iostream>
#include<iomanip> // 带上这个库，我们就可以用 setprecision 指定输出的精度
using namespace std;

int main(){

	double a = 10000.1, b = 10000;
	double c = a - b;
	cout<<" relative error = "<<setprecision(15)<< c/0.1 -1<<endl; // c 的理论值为 0.1， 所以 c/0.1 -1 就是相对误差
	// setprecision(15)指定输出 15 位有效数字

	return 0;
}
\end{lstlisting}
运行以后，会得到
\begin{lstlisting}
	relative error = 3.63797880709171e-12,
\end{lstlisting}
注意，$e-12$ 是指 $10^{-12}$。
根据我们上面的理论分析，相对误差的上限约为
\begin{equation}
	\frac{ 10000 }{0.1} 10^{-16} = 10^{-11},
\end{equation}
所以结果与我们的理论预测可以认为是一致的。
你可以试着改动这个程序，将 a,b 中的 0 多加几个，观察程序输出的相对误差，体会体会我们的分析结果的力量。

\section{第1部分作业}
1. 在　c++　中，求余运算符号为　\%，例子：
\begin{lstlisting}
    a = 3%2;
\end{lstlisting}
即用 3　除以　2　得到的余数，即　１　赋值　ａ。
请用我们讲的知识，编写程序，完成以下任务
\begin{itemize}
\item 要求用户输入两个正整数 a,b
\item 在屏幕上输出　a 除以　b 的余数。
\end{itemize}

\vspace{1cm}

*2. 质数的定义为：只能被1和自己整除的正整数。编写程序输出100以内的所有质数。

\vspace{1cm}


*3. 编写程序，完成如下任务：\\
\begin{itemize}
\item 要求用户输入 7 个浮点数
\item 自动从大到小排序，并输出到屏幕
\end{itemize}

**4. $e^x$　可以做泰勒展开
\begin{equation}
	e^x = \sum^\infty_{n=0} \frac{x^n}{n!},
\end{equation}
可以指定当　$ |x^N / N!| < \epsilon$　（$\epsilon$　是一个很小的数，比如，可以取　$\epsilon = 10^{-6}$）时，将前面　$N$　项的和作为　$e^x$　的近似：
\begin{equation}
	e^x \approx \sum^{N-1}_{n=0} \frac{x^n}{n!},
\end{equation} 

也可以将求和顺序反过来
\begin{equation}
	e^x \approx \sum^{0}_{n=N-1} \frac{x^n}{n!},
\end{equation} 
请问哪一种算法可以更有效地避免舍入误差的影响，得到更好的精度？
试编程，分别用这两种算法计算　$e$　的值。\\
\vspace{1cm}
提示：c++　中如果包括　cmath 包，即
\begin{lstlisting}
    #include<cmath>
\end{lstlisting}
即可在代码中调用　e 指数函数，即　exp(１)　表示　$e^１$，作为精确值，与上面的两种算法的计算结果相比较。

*5. 请参考百度词条“利萨如图形”,并用gnuplot做出任意一个利萨如图形，并输出到一个eps文件。（要求提交gnuplot脚本，并说明你的利萨如图形的方程。提示：gnuplot脚本中，＃后面的内容表示注释。）

**6. 如果我们需要做求和，
\begin{equation}
s_{up} = \sum^N_{n=1} 1/n,
\end{equation}
我们可以按上式的顺序求和，也可以用逆序求和
\begin{equation}
s_{down} = \sum^1_{n=N} 1/n,
\end{equation}
１）请问哪一种方式更精确？\\
２）编写程序，计算 $log_{10} | (s_{up} (N) - s_{down} (N) )/ s_{down} (N) |$，$N=10,100,1000,...,10^{10}$。
将这10组数据输出到文件。
３）用上面的数据，在gnuplot中画出折线图，用 terminal ps，输出到一个ps文件。